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SYSTEMS AND METHODS FOR CORRECTING 
INHOMOGENEITY IN IMAGES 

BACKGROUND OF THE INVENTION 

[0001] The invention disclosed and claimed herein is generally 
directed to a method for correcting spatial inhomogeneity or nonuniformity of spatial 
intensity in an acquired magnetic resonance (MR) or other medical diagnostic image. 
More particularly, the invention is directed to a correction method of such type 
wherein the primary component of inhomogeneity is slowly varying. 

[0002] In many areas of imaging including MR and computed 
tomography (CT), acquired images are corrupted by slowly varying multiplicative 
inhomogeneities or nonuniformities in spatial intensity. Such nonuniformities can 
hinder visualization of the entire image at a given time, and can also hinder automated 
image analysis. Such inhomogeneity is a particular concern in MR when single or 
multiple surface coils are used to acquire imaging data. The acquired images 
generally contain intensity variations resulting from the inhomogeneous sensitivity 
profiles of the surface coil or coils. In general, tissue next to the surface coil appears 
much brighter than tissue far from the coil. Spatial intensity variations introduced by 
surface coil nonuniformity hinders visualization because one cannot find a 
window/level adjustment to encompass the entire field of view. When such images 
are filmed, the operator tries to select a setting which covers most of the features of 
interest. Furthermore, uncorrected image inhomogeneity makes it difficult to perform 
image segmentation and other aspects of image analysis. 

[0003] An example of the problem is spine imaging, wherein one or 
more surface coils are placed behind a patient. If the central spinal canal is filmed 
optimally, tissue structure behind the vertebral column may be overamplified and may 
become so bright that no tissue detail can be seen. At the same time, tissue in front of 
the vertebral column may be so dark that image detail in that area is also obscured. 
Therefore, in order to optimally display and film the entire image, the signal variation 
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due to the inhomogeneous sensitivity profile of the surface coil needs to be corrected. 
Surface coil image signal intensities generally represent the product of (1) precessing 
magnetization of the body tissue or other object being imaged, and (2) the sensitivity 
profile of the surface coil. Accordingly, various intensity correction algorithms have 
been devised to correct surface coil images by dividing out an estimate of the surface 
coil's sensitivity profile. Thus, if the observed or acquired MR image signal is defined 
in a spatial domain for a pixel location (x,y) by a function g(x,y) then 
g(x,y)=h(x,y)*f(x,y) + n(x,y), where * represents multiplication, h function., f function 
and n represent the coil profile function, a corrected function, and the imaging noise, 
respectively. More specifically, the corrected function is a function defining an image 
which is substantially free of distortion resulting from the inhomogeneity. Thus, the 
problem is to determine both h and f, given only the measured or acquired function g 
in the presence of n. However, if the function h can be determined which reasonably 
represents the inhomogeneity distortion, then the function f can be readily computed 
from f = [g*h]/[h*h+yi], which is known in the art as Werner filter solution, where \\f\ 
is a regularization parameter corresponding to the reciprocal of signal to noise ratio. 

[0004] The distortion arising from use of surface coils generally 
varies slowly over space. An important class of prior art solutions to the above 
problem is based on this assumption. In accordance therewith, a low pass filtering 
operation is applied to g. The resulting function, represented as LPF[g], does not 
contain high frequency components and is taken as an estimate of distortion function 
h. An estimate of the function f is then obtained by dividing g by LPF[g], i.e., f = 
g/LPF[g]. However, for this class of methods to be effective, g must not contain sharp 
intensity transitions. Unfortunately, in MR imaging an air-lipid interface usually 
contains sharp intensity transitions which violate the basic assumption made in the 
method, i.e., that the low frequency content in the scene being imaged is solely due to 
h. Significant air-lipid interferences will generally be encountered, for example, at the 
edges of an organ, i.e., at the boundary between the organ and an air-space or cavity. 

[0005] To overcome the above deficiency in low pass filtering 
correction at the edge or boundary of an organ or other tissue structure, certain hybrid 
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filtering techniques have been developed. Some of such techniques are set forth in the 
following references: Surface Coil MR Imaging of the Human Brain with an Analytic 
Reception Profile Correction, JMRI 5, 139-144, by S. E. Moyher, D. B. Vigeron, and 
S. J. Nelson; Phased Array Detectors and an Automated Intensity Correction 
Algorithm for High Resolution MR Imaging of the Human Brain, JMRI (1995), by L 
L. Wald, L. Carvajal, S. E. Moyher, S. J. Nelson, P. E. Grant, A. J. Barkovich, and D. 
B. Vingeron; Phased Array Image Intensity Correction: An Algorithm to Remove 
Intensity Variations in MR Images Resulting from the Inhomogeneous Sensitivity 
Profiles of Phased Array Surface Coils, a Master's thesis by J. Murakami (1995), 
University of Washington, Seattle, and in U.S. Patent 5,943,433. 

[0006] Intensity inhomogeneities can influence visualization of 
digital images. If there are multiple kinds of inhomogeneities in a digital image, 
several difficulties arise while compensating for such inhomogeneities. As an 
illustration, if there are two different kinds of inhomogeneities in a digital image, a 
designed algorithm may not be able to properly multiple inhomogeneities since the 
algorithm is tuned to mitigate one kind of inhomogeneity. Inadequate compensation 
can result in unsatisfactory visual appearance of the compensated image. If any 
specific scale is used for intensity correction, it is inadequate due to the reasons stated 
above. On the other hand, if an entire image is used for the intensity correction, it 
results in a flat image without any contrast. Moreover, if a normalization based on 
average intensity is performed on digital images, inconsistent results can occur 
between different slices of the digital images. 

BRIEF DESCRIPTION OF THE INVENTION 

[0007] In one aspect a method for generating an estimate of 
inhomogeneity is provided. The method includes generating a first estimate of 
inhomogeneity, generating a second estimate of inhomogeneity, and generating a final 
estimate of inhomogeneity using at least the first and second estimates. 

[0008] In another aspect, a magnetic resonance imaging (MRI) 
system includes a main magnet configured to generate a substantially uniform 
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magnetic field, a radio frequency pulse generator configured to excite the magnetic 
field, a gradient field generator configured to generate gradients extending in different 
directions in the magnetic field, a receiver configured to receive magnetic field 
magnetic resonance (MR) signals representative of an object, and a computer 
operationally coupled to the receiver. The computer is configured to generate a first 
estimate of inhomogeneity, generate a second estimate of inhomogeneity, and generate 
a final estimate of inhomogeneity using at least the first and second estimates. 

[0009] In another aspect, a computer readable medium is encoded 
with a program configured to instruct a computer to generate a first estimate of 
inhomogeneity, generate a second estimate of inhomogeneity, and generate a final 
estimate of inhomogeneity using at least the first and second estimates. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0010] Figure 1 is a simplified block diagram of a magnetic 
resonance imaging (MRI) system which incorporates systems and methods for 
correcting inhomogeneity in images. 

[0011] Figure 2 is a flowchart of an embodiment of a method for 
correcting inhomogeneity in images. 

[0012] Figure 3 shows images of a brain of a subject scanned using 
the MRI system. 

DETAILED DESCRIPTION OF THE INVENTION 

[0013] The herein described systems and methods provide 
equalization in images having more than one kind of inhomogeneities. A final 
homogeneity estimate or digital values h(x,y) incorporate a combination of a specific 
level of smoothness and an entire level of smoothness. For instance, a first 
inhomogeneity estimate hi(x,y) and a second inhomogeneity estimate h 2 (x,y) are 
combined into the final inhomogeneity estimate h. This combination can be 
accomplished using a suitable blending of the specific scale image and an acquired 
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image I a . Moreover, the herein described systems and methods provide improvements 
in the normalization to compensate for potential variations in intensities of digital 
images. Additionally the herein described systems and methods automatically 
determine a threshold for compensating smoothing of strong image features such as 
those appearing at tissue-air interfaces. 

[0014] Figure 1 is a simplified block diagram of a magnetic 
resonance imaging (MRI) system 10 which incorporates systems and methods for 
correcting inhomogeneity in images. MRI system 10 is illustrated as including a 
scanner 12 coupled to circuitry for acquiring and processing discrete pixel data. 
Scanner 12 includes a support structure 14 in which a physical subject 16, such as a 
human being, may be placed for acquiring images representative of internal features, 
such as tissues, fluids and so forth. Scanner 12 includes an electromagnet 
arrangement 18 for producing an electromagnetic field. Excitation and sensing coils 
20 are provided within scanner 12 for exciting gyromagnetic materials within subject 
16 and for sensing emissions from the materials. 

[0015] Signals sensed by coils 20 are encoded to provide digital 
values representative of the excitation signals emitted at specific locations within the 
subject, and are transmitted to signal acquisition circuitry 22. Signal acquisition 
circuitry 22 also provides control signals for configuration and coordination of fields 
emitted by coils 20 during specific image acquisition sequences. Signal acquisition 
circuitry 22 transmits the encoded image signals to a signal processing circuit 24. 
Signal processing circuit 24 executes pre-established control logic routines stored 
within a memory circuit 26 to filter and condition the signals received from signal 
acquisition circuitry 22 to provide 50 (see figure 2) digital values ga(x,y), such as, for 
instance, intensity values, representative of each pixel at location (x,y) in the acquired 
image I a . The digital values g a (x,y) are then stored in memory circuit 26. Signal 
processing circuitry 22 executes the methods for correcting inhomogeneity in the 
acquired image I a by using the digital values g a (x,y). Although I a refers to an entire 
image, and g a (x,y) refers to a particular pixel, I a and g a (x,y) are sometimes used 
interchangeably herein. Additionally, to ease readability, ga is also used to represent 
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g a (x ? y). This convention is also maintained in reference to other images and 
respective pixel values herein. 

[0016] Signal processing circuit 24 receives configuration and 
control commands from an input device 28 via an input interface circuit 30. Input 
device 28 will typically include an operators station and keyboard for selectively 
inputting configuration parameters and for commanding specific image acquisition 
sequences. Signal processing circuit 24 is also coupled to an output device 32 via an 
output interface circuit 34. Output device 32 will typically include a monitor or 
printer for generating reconstituted images based upon the image enhancement 
processing carried out by signal processing circuit 24. 

[0017] It is noted that, while in the present discussion reference is 
made to discrete pixel images generated by an MRI system, the systems and methods •* 
described herein are not limited to any particular imaging modality. Accordingly, the 
systems and methods may also be applied to image data acquired by X-ray systems, 
digital radiography systems, positron emission tomography (PET) systems, and 
computed tomography (CT) systems, among others. It should also be noted that in the 
embodiment described, signal processing circuit 24, memory circuit 26, and input and 
output interface circuits 30 and 34 are included in a programmed digital computer 36. 
However, circuitry for carrying out the techniques described herein may be configured 
as appropriate coding in application-specific microprocessors, analog circuitry, or a 
combination of digital and analog circuitry. Computer 36 includes a device 38, for 
example, a floppy disk drive, CD-ROM drive, DVD drive, magnetic optical disk 
(MOD) device, or any other digital device including a network connecting device such 
as an Ethernet device for reading instructions and/or data from a computer-readable 
medium 40, such as a floppy disk, a CD-ROM, a DVD or an other digital source such 
as a network or the Internet, as well as yet to be developed digital means. In another 
embodiment, computer 36 executes instructions stored in firmware (not shown). 
Computer 36 is programmed to perform functions described herein, and as used 
herein, the term computer is not limited to just those integrated circuits referred to in 
the art as computers, but broadly refers to computers, processors, microcontrollers, 
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microcomputers, programmable logic controllers, application specific integrated 
circuits, and other programmable circuits, and these terms are used interchangeably 
herein. 

[0018] Figure 2 is a flowchart of an embodiment of a method for 
correcting inhomogeneity in images. The method includes distinguishing between 
structure and non-structure in the acquired image I a by performing 52 a mask 
operation. The mask operation produces a mask image I maS k having digital values 
gmask(x,y). In the mask operation, digital value of a pixel is set to 1 if the 
corresponding pixel (x,y) in image I a is structure and is set to 0 if the corresponding 
pixel (x,y) in image I a is non-structure. 

[0019] The method distinguishes structure from non- structure as 
follows. A normalized image I n0 rm is formed from the acquired image I a . The 
acquired image I a is normalized by reading the digital values g a (x,y) and scaling the 
values over a desired dynamic range. A gradient image Igr a d is computed from a 
blurred or smoothed version of the normalized image I nor m. An example of a 
smoothing operation includes boxcar smoothing. In boxcar smoothing, a boxcar filter 
smoothes the normalized image I nonT1 by computing the average of a given 
neighborhood. The kernel is separable and its length, such as, for instance, 3 pixels, is 
variable. The kernel is moved horizontally and vertically along the normalized image 
Inorm until each pixel of the image has been processed. The smoothing operation 
results in an image I s having digital values gs(x,y). 

[0020] For example, at every pixel of the image I s , computed are an 
X gradient component, a Y gradient component, and a resultant gradient as the 
maximum of the two components. If desired, the direction of the gradient using 
arctan(Y gradient component/X gradient component) can be computed. While several 
techniques may be employed for this purpose, in one embodiment, 3X3 Sobel 
modules or operators are employed. One module is used for identifying the X 
gradient component, while another module is used for identifying the Y gradient 
component of each pixel of the image Is. In this process, the modules are 
superimposed over an individual pixel of interest in the image I s , with the pixel of 
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interest situated at the central position of the 3X3 module. The digital values g s (x,y) 
located at element locations within each module are multiplied by the scalar value 
contained in the corresponding element, and the resulting values are summed to arrive 
at the corresponding X and Y gradient components. 

[0021] With these gradient components thus computed, the gradient 
magnitude G mag and gradient direction Gdi r are computed. The gradient magnitude 
G mag for each pixel is equal to the higher of the absolute values of the X and Y 
gradient components for the respective pixel. The gradient direction is determined by 
finding the arctangent of the Y gradient component divided by the X gradient 
component. For pixels having an X gradient component equal to zero, the gradient 
direction is assigned a value of n/2. The values of the gradient magnitudes Gm ag and 
gradient directions G d i r for each pixel are saved in memory circuit 26. 

[0022] It should be noted that alternative techniques may vbe 
employed for identifying the X and Y gradient components and for computing the 
gradient magnitudes and directions. For example, in place of Sobel gradient modules, 
other modules such as the Roberts or Prewitt operators are employed. Moreover, the 
gradient magnitude G ma g may be assigned in other manners, such as a value equal to 
the sum of the absolute values of the X and Y gradient components or a value equal to 
a square root of a sum of squares of the X and Y gradient components. 

[0023] A gradient histogram is created specifically based upon the 
gradient magnitudes G mag . The gradient histogram is a bar plot of specific populations 
of pixels having specific gradient magnitudes G mag . These gradient magnitudes G mag 
are indicated by positions along a horizontal axis, while counts of the pixel 
populations for each gradient magnitude G mag are indicated along a vertical axis, with 
each count falling at a discrete level. The resulting bar graph forms a step-wise 
gradient distribution curve. It is noted that in an actual implementation, the gradient 
histogram need not be represented graphically, but may be functionally determined by 
signal processing circuit 24 operating in cooperation with values stored in memory 
circuit 26. 
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[0024] The gradient histogram is used to identify an initial gradient 
threshold (IGT) for separating structure from non-structure. The IGT is set at a 
desired gradient magnitude level G mag . Pixels having gradient magnitudes G mag at or 
above the IGT are considered to meet a first criterion for defining structure, while 
pixels having gradient magnitudes G mag lower than the IGT are initially considered 
non-structure. The IGT used to separate structure from non- structure is set by an 
automatic processing or "autofocus" routine. However, alternatively, IGT is set by 
operator intervention, via input interface circuit 30 or a value identified through the 
automatic process described above may be overridden by the operator. 

[0025] The process for identification of the IGT is as follows. The 
IGT is conveniently set to a value corresponding to a percentile of the pixel 
population, such as 30 th percentile. The location of the IGT along the horizontal axis 
is thus determined by adding pixel population counts from the left-hand edge of the 
gradient histogram, adjacent to the vertical axis and moving toward the right. Once a 
desired percentile value is reached, the corresponding gradient magnitude G mag is the 
value assigned to the IGT. It is noted that in an alternative embodiment, the 
smoothing and the normalization operations that are described above are not 
performed. 

[0026] The mask image I mask is isotropically smoothed in the non- 
structure region of the image and anisotropically smoothed in the structure region of 
the image to obtain 54 a filtered image If with digital values go(x,y). Boxcar 
smoothing, which is an example of isotropic smoothing, is described above. In 
anisotropic smoothing, the structure region is filtered to extract domain orientation 
information. Anisotropic smoothing involves iteratively filtering the structure by a 
smoothing kernel, such as a 3X1 smoothing kernel, along a dominant direction in a 
given neighborhood of the image I mask , which would be the direction of majority of the 
local minimum variances in that neighborhood. The iterations are performed a set 
number of times, such as, for instance, 3 times. 

[0027] Each iteration is accomplished using the following steps. The 
structure region is scanned and a local orientation map is obtained by assigning one of 
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four orientation numbers, for instance, assigning 1 for 45 degrees, 2 for 135 degrees, 3 
for 90 degrees, and 4 for 0 degrees. The structure region is scanned again and the 
dominant orientation at any point is determined by counting the number of different 
orientations in a neighborhood of the image I mask . An orientation getting the 
maximum counts is the dominant orientation. 

[0028] As a further refinement, both the dominant orientation and its 
orthogonal orientation are used to make a consistency decision. The consistency 
decision is made if one of the following conditions is met: (1) An orientation getting 
maximum counts is greater than a pre-specified percentage, for instance, 0.67 percent, 
of the total counts in a neighborhood of the image I mas k and its orthogonal orientation 
gets the minimum counts, (2) An orientation getting maximum counts is greater than a 
smaller percentage, for example, 0.44 percent, of the total counts in a neighborhood of 
the image Imask, its orthogonal orientation gets the minimum count, and the ratio of the 
dominant count and its orthogonal count is greater than a pre-specified number, for 
example, 5, and (3) The ratio of a dominant count to its orthogonal count is greater 
than 10. Smoothing is performed along a dominant orientation. 

[0029] The acquired image I a and the filtered image If are blended 56 
using different parameters corresponding to the structure and the non-structure regions 
to obtain a blended image lb with digital values g(x,y). As an example, the acquired 
image I a and the filtered image If are blended using an equation: 

g(x,y) - a(go(x,y) - g a (x,y)) + g a (x,y), if g mask (x,y) = 0 
else g(x,y) = p(g G (x,y) - g a (x,y)) + g a (x,y), if gmask(x,y) = 1, 

...(1) 

where a is an operator selected parameter such that 0 < a<l, and where 
p is an operator selected parameter such that 0 < p<l. 

[0030] The method for correcting inhomogeneity in images further 
includes calculating 58 a threshold value T from the digital values g 0 (x,y) of the 
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filtered image If. Alternatively, the threshold value T is calculated directly from the 
mask image I mas k- From the filtered image I f , a foreground intensity (FI) is computed 
as an average of the structure regions and a background intensity (BI) is computed as 
an average of the non-structure regions having intensity less than FI. From the non- 
structure regions with intensities less than FI, a standard deviation (SD) is also 
calculated. The threshold value T is obtained from FI, BI, and SD. As an example, 



If (SD/BI) < 0.2, T = FI * (BI/FI +0.2); 
Else if (SD/BI) < 0.5, T = FI * (BI/FI +0.1); 
Else if (SD/BI) < 1.0, T = BI; 

Else if ((SD/BI) > 1.0) AND (BI/FI) < 0.2, T = FI * 0.2; 
Else { 

T = FI*(BI/FI^0.1); 
IfT<0.0,T = FI*0.1; 



It is noted that the values 0.1, 0.2, and 0.5 in equation (2) are 
exemplary and that any value between 0 and 1 can be used instead. 



[0031] The method includes determining 60 an inhomogeneity image 



linhomogendty having digital values h(x,y) from the digital values g(x,y) of the blended 
image lb. The inhomogeneity image Iinhomogeneity corrupts a desired image Idesired having 
digital values f(x,y). The relationship between the digital values h(x,y), g(x,y), f(x,y), 
and additive noise n is provided by 



[0032] The inhomogeneity image linhomogendty is obtained as follows. 
The digital values g(x,y) are shrunk by a factor equal to or greater than 1. As an 
example, the shrink function uses average digital values in non-overlapping integer 



neighborhoods of a particular size, such as, 3X3, to create a shrunk image I S hmnk with 
digital values g S hrunk(x,y). It is noted that other types of shrink functions can be used. 



} 



•••(2) 



g = hf+n 



-(3) 
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A thresholding function is performed on the digital values g S hrunk(x,y) based on the 
threshold value T. For example, if the digital value g s hrunk(x,y) is greater than the 
threshold value T, g S hrunk(x,y) is set to 1 and if g s hrunk(x,y) is less than T, g S hrunk(x,y) is 
assigned a value of 0. The thresholding function results in digital values 
Threshold(g shrunk (x,y)). The digital values g shmnk (x,y) and Threshold(g shr unk(x,y)) are 
filtered by using, for example, a low pass filter (LPF), to provide digital values, 
LPF(gshrunk(x,y)) and LPF[Threshold(gs h runk(x,y))] respectively. LPF removes higher 
frequency components associated with sharp intensity transitions. LPF commences by 
taking a transform, such as a fast Fourier transform (FFT) or a Discrete Cosine 
Transform (DCT), of gshrunk(x,y) and Threshold(g shr unk(x,y). Respective transform 
components are then multiplied by coefficients predetermined in accordance with a 
filter operation, such as a Gaussian filter operation. Inverse functions, such as an 
inverse FFT, are then computed for respective coefficient multiplication's, to 
determine LPF(g shrun k(x,y)) and LPF[Threshold(g shrunk (x,y))] . 

[0033] Digital values h s (x,y) are obtained as follows. 

h s (x,y) = LPF(g shrunk (x,y))/LPF[Threshold(g shrunk (x,y))] 

...(4) 

First inhomogeneity estimates or digital values hi(x,y) are obtained 
from the digital values h s (x,y) by expanding h s (x,y). The digital values h s (x,y) are 
expanded by implementing a process such as bilinear interpolation and for 3- 
dimensional digital values, trilinear interpolation is used to expand the 3-dimensional 
digital values. Such a process is described, for example, in Numerical Recipes in C, 
The Art of Scientific Computing, by William H. Press, Brian P. Flannery, Saul A. 
Tenkolsy, and William T. Vetterling; Cambridge University Press (1988). 

[0034] Second inhomogeneity estimates or digital values h 2 (x,y) are 
obtained by dividing the digital values g(x,y) by Threshold(g(x,y)). Alternatively, the 
digital values h 2 (x,y) are obtained by dividing the digital values LPF(g(x,y)) by 
LPF[Threshold(g(x,y))]. In yet another alternative embodiment, the digital values 
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h 2 (x,y) are obtained by dividing the digital values g 0 (x,y) by Threshold(g 0 (x,y)). In 
still another alternative embodiment, the digital values h 2 (x,y) are obtained by 
dividing the digital values LPF(go(x,y)) by LPF[Threshold(go(x,y))}. In another 
embodiment, the digital values h 2 (x,y) are obtained by smoothing the mask image 
Imask at a different level, by obtaining the blended image I b with the digital values 
g(x,y), and by dividing g(x,y) by Threshold(g(x,y)). In yet another embodiment, the 
digital values h 2 (x,y) are obtained by smoothing the mask image I mask at a different 
level, by obtaining the blended image lb with the digital values g(x,y), and by dividing 
LPF(g(x,y)) by LPF[Threshold(g(x,y))]. In stijl another embodiment, the digital 
values h 2 (x,y) are obtained without filtering, such as low pass filtering, g(x,y) or 
go(x,y). 

[0035] Final inhomogeneity estimates or the digital values h(x,y) of 
the inhomogeneity image Ihompgeneity are obtained as follows. 

h(x,y) = h 1 + (112-1100, ...(5) 

where O<0<1. 

In an alternative embodiment, 

h(x,y) = Q { h t + 8 2 h 2 + . . . + 8 N h N , . . .(6) 

where 0<9i, 9 2 , . . 0 N <1 and 0 r + 0 2 . + . . . + 0 N = 1 

In the alternative embodiment, h„ where i = 1 ... .N, is obtained by using 
different filters, such as, lowpass, highpass, bandpass, and Gaussian filters, other 
functions, such as, a shrink function, an average function, a median function, weighted 
polynomials, a minima function, a maxima function, and any combination of the 
filters and the functions. For example, 

hj = [minima of digital values within a neighborhood of 
g(x,y)]/[Threshold(minima of digital values within the neighborhood of g(x,y))]. 

...(7) 
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As another example, 

hj =5= [average of digital values within a neighborhood of 
g(x,y)]/[minima of digital values within the neighborhood of g(x,y)]. 

•••(8) ' 

As yet another example, 

hj = Gaussian filter(g(x,y))/(minima of digital values within a 
neighborhood of g(x,y)). . . . (9) 

As another example, each h, is obtained using a different LPF. For 
instance, hj is obtained using an LPF filtering a first frequency range and h 2 is 
obtained using an LPF filtering a second frequency range. It is noted that to compute 
hi, any of or any combination of go, g a , gshrunk and g can be used. It is also noted that 
to compute the threshold value T, any of or any combination of g 0 , g a , gshrunk and g can 
be used. Therefore as used herein g m refers to any of go, g a , gshrunk and g. 

[0036] The digital values f(x,y) of the desired image I<j e sired are 
obtained 62 as follows. 

[0037] f=[gh]/[hh + «, ...(10) 

where § > 0, such as, for instance, 1 . 

It is noted that by choosing a large 0, such as, for example, 0.3, by 
shrinking hj by a large amount, and by performing the anisotropic and isotropic 
smoothing at a high level, high computational efficiency is achieved. 

[0038] The method includes determining 64 gradient magnitudes of 
the digital values g(x,y) of the blended image lb and determining 66 gradient 
magnitudes of the digital values f(x,y) of the desired image Idesired. A process of 
determining gradient magnitudes is described above. Digital values fjv of a 
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normalized image I N are obtained by computing 68 a normalization factor and 
normalizing 70 the desired image Idesired as follows. 

f N = [{f(grad(g))}/grad(f)]M, ...(11) 

where M is a multiplier greater than or equal to 1 , such as, for example, 
10, [(grad(g))/grad(f)]M is the normalization factor, and grad is a gradient operation 
that determines gradient magnitudes as described above. 

[0039] The method includes restricting a ratio f N /g. The ratio is 
restricted as follows. 

If(f N /g)>L K ,fN = gLK, 

where L K is an empirically determined factor whose value is greater 
than or equal to 1, such as, for example, 6. The restriction usually results in an image 
that does not have large areas of dark regions and large areas of bright regions by 
compensating for a high level of isotropic and anisotropic smoothing. It is noted that 
although the methods and systems described herein use 2-dimensional digital values, 
the methods and systems can use 3-dimensional digital values instead. 

[0040] Figure 3 shows images of a brain of subject 16. Image 90 is 
an example of the blended image lb having the digital values g(x,y) and image 92 is an 
example of the image Idesired, having the digital values f(x,y), obtained after executing 
the methods for correcting inhomogeneity in images. Image 92 shows lesser 
inhomogeneities than inhomogeneities present in image 90 and so is an improvement 
over image 90. 

[0041] While the invention has been described in terms of various 
specific embodiments, those skilled in the art will recognize that the invention can be 
practiced with modification within the spirit and scope of the claims. 
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